An Efficient and Accurate Car-Parrinello-like 
Approach to Born-Oppenheimer Molecular Dynamics 
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We present a new method which combines Car-Parrinello and Born-Oppenheimer molecular dy- 
namics in order to accelerate density functional theory based ab-initio simulations. Depending on 
the system a gain in efficiency of one to two orders of magnitude has been observed, which allows 
ab-initio molecular dynamics of much larger time and length scales than previously thought feasible. 
It will be demonstrated that the dynamics is correctly reproduced and that high accuracy can be 
maintained throughout for systems ranging from insulators to semiconductors and even to metals 
in condensed phases. This development considerably extends the scope of ab-initio simulations. 
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Density functional theory (DFT) [l[ based ab-initio 
molecular dynamics in which interactions are com- 
puted on the fly from electronic structure calculations, 
has been very successful in describing a large variety of 
physical phenomena and has proven its relevance in many 
fields. However, its computational cost has limited the 
attainable length and time scales in spite of substantial 
progress. For a while it was believed that the devel- 
opment of linear scaling methods 0, 0, HI could have 
offered a solution. Unfortunately, the crossover point 
at which linear scaling methods become advantageous 
has remained fairly large, especially if high accuracy is 
needed. 

On the other hand, it would be very desirable to accel- 
erate ab-initio simulations with up to thousands of atoms, 
such that simulations as long as tens or even hundreds 
of picoseconds can be routinely performed, thus making 
completely new phenomena accessible to ab-initio simula- 
tions. Born-Oppenheimer molecular dynamics (BOMD) 
simulations, in which at each time step the DFT func- 
tional is fully minimized, do not seem to offer much room 
for further improvement. For this reason another direc- 
tion has recently been followed to improve the efficiency 
at current system sizes. In the spirit of Car-Parrinello 
molecular dynamics (CPMD) [2], some form of dynam- 
ics for the electronic degrees of freedom is implemented, 
which automatically maintains the system close to the 
BO surface, but at variance with the original proposal in 
a localized orbital representation 0, 0, H@] ■ The acceler- 
ation stems from the ability to reduce or fully bypass the 
self-consistency cycle. However, just like in CPMD, these 
methods suffer from rather short integration time steps. 
In this Letter we present a new method which overcomes 
this limitation and combines the accuracy and long time 
steps of BOMD with the efficiency of CPMD. 

We consider the general case in which the DFT Kohn- 
Sham orbitals are expanded in a non-orthogonal basis. 
Let M be the dimension of the Hilbcrt space and S the 
M x M overlap matrix. As is customary in quantum 



chemistry we arrange the expansion coefficients of the N 
lowest occupied orbitals in a rectangular M x N matrix 
C. The density matrix is then written as P = CC T 
and obeys the idempotency condition P = PSP. Main- 
taining the idempotency of P is crucial in any electronic 
structure calculation and is one of its algorithmic chal- 
lenges. The potential energy surface on which the ions 
move is defined by the minimum of an appropriately cho- 
sen energy functional -Edft [C, Rj], which we express as 
a functional of C and a function of the ionic coordinates 
R/. In this notation the Born-Oppenheimer equation of 
motions reads as follows: 



M/Rj = - V/ nun £ D ft [C, Rj] , 



(1) 



where the search for the minimum is restricted to the 
C's that satisfy the orthonormality condition C T SC = I, 
which is equivalent to imposing the idempotency condi- 
tion on P. The forces of Eq. (Q} can be divided into three 
contributions, the Hellmann-Feynman forces 0, [ll|, the 
Pulay forces [l2j, which are present whenever the basis 
set depends on the ionic positions, and a residual term 
[l3| which is non-zero except when full self-consistency 
is reached. This term leads to poor energy conservation 
in BOMD unless a very tight convergence criterion is im- 
posed. In Car-Parrinello-like approaches this is circum- 
vented by the design of a coupled electron-ion dynamics 
which maintains the system close to the BO surface, but 
at the cost of small time steps. 

Based on ideas of the original Car-Parrinello approach 
we design an improved dynamics for the coupled system 
of electrons and ions. Contrary to the Car-Parrinello 
scheme we will not write an explicit equation of motion 
for the C's, but rather an integration scheme for the elec- 
tronic degrees of freedom. The knowledge of the previ- 
ous K values of C (£„_/), where I € [1,K], determines 
the value of C (t n ), such that at any instant of time the 
C's are as close as possible to the instantaneous ground 
state. As for the short-term integration of the electronic 
degrees of freedom, accuracy is crucial, a highly accurate 
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and efficient algorithm is required. Here we have cho- 
sen the always stable predictor-corrector (ASPC) method 
14j of Kolafa. This method was originally designed to 
deal with classical polarization, so that additional care 
must be taken that during the evolution the idempotency 
condition is satisfied. The predictor 
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uses the extrapolated contra-covariant density matrix PS 
[l5| as an approximate projector on to the occupied sub- 
space C(<„_i). In this way, we take advantage of the 
fact that the physically relevant contra-covariant density 
matrix PS evolves more smoothly and is therefore easier 
to predict than C. This is followed by a corrector step 
to minimize the error and to further reduce the deviation 
from the ground state. The corrector 



C (t n ) = w MIN [C p (*„)] + (1 - uj)C p (t n ) . 
K 

with lo 
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(3) 



consists of a single minimization step MIN [C p (t„)] of a 
properly selected minimization procedure. Alternatively, 
the predictor can also be repeatedly applied. In this case 
the ground state is even more closely approached, but 
in general this is not necessary. The numerical coeffi- 
cients of Eq. © were selected by Kolafa in order to en- 
sure time-reversibility up to 0(h K+2 ), while co was cho- 
sen to guarantee a stable relaxation towards the mini- 
mum. Because the energy is invariant under a unitary 
transformation within the subspace of occupied orbitals 
C, it must be ensured that this gauge transformation is 
not strongly changed by MIN [C p (tn)], as in this case 
continuity between the C's may be lost. Moreover the 
minimization scheme must be very efficient in bringing 
the system close to the ground state and preserving the 
idempotency of the density matrix. For these reasons we 
have chosen the orbital transformation (OT) method of 
VandeVondele and Hutter [ijj]. Inspired by the form of 
the exponential transformation 17[ an auxiliary variable 
X is introduced, to parameterize the occupied orbitals 

C (X) = C p (t n ) cos (U) + XLT 1 sin (U) (4a) 

1 /2 

where U = (X T SX) and the variable X has to obey 
the linear constraint: 



X T SC P (t n ) = 



(4b) 



Under this condition C (X) leads to an idempotent den- 
sity matrix for any choice of X, provided that the ref- 
erence orbitals C p (t n ) are orthonormal. Thus a finite 
minimization step along the gradient direction will still 



fulfill the idempotency constraint exactly. Due to the 
linear constraint the minimization with respect to X is 
performed in an auxiliary tangent space. Because this is 
a linear one no curved geodesies must be followed, as is 
the case for variables such as C, which are nonlinear ly 
constrained. In this way large minimization steps can 
be taken, especially if a good preconditioner is used. In 
fact using a very efficient, idempotency conserving direct 
minimizer like OT has been decisive for the success of this 
approach. Since the ASPC integrator preserves the or- 
thonormality constraint only approximately, it occasion- 
ally has to be explicitly enforced, either by a Cholesky 
decomposition or by a few purification iterations [l8j |. 

Having obtained the new wavefunction we can now 
evaluate the energy and the nuclear forces, which are de- 
rived from the following approximate energy functional: 



Epc\p] = Tr [C T H[p'P] C]-±JdrJ dr 
- Jdr V xc [p p ]p p + E xc [p p ] + E H , 



, p p (r)p p (r>) 



(5) 



where p p (v) is the density associated with C p (t n ). Epc[p] 
can be thought of as an approximation to the Harris- 
Foulkcs functional [191 ] and maintains the predictor- 
corrector flavor of the present work. The validity of 
Epc[p] depends only on the efficiency of the minimizer 
and on the quality of the propagation scheme. The ionic 
forces are calculated by evaluating the analytic gradient 
of Epc[p], with respect to the nuclear coordinates. How- 
ever, as Ap = p — p p ^ 0, besides the usual Hellmann- 
Feynman and Pulay forces an extra term appears: 



dr 



dVxc[p p ] 
dpP 



Ap + Vn[Ap] 



(6) 



where p is the corrected density evaluated using C(i„) 
and p p is the predicted density calculated from C p (t n )- 
However, as a single minimization step is performed, 
C(t n ) is only an approximate cigenfunction of H[p p ] 
within the subspace spanned by the finite basis set used 
[20| . This leads to an insignificant error in the forces, 
provided that C(i„) is very close to the ground-state. 

The ability of this dynamics to maintain the system on 
the BO surface may vary considerably. It is almost per- 
fect in systems like water, but somewhat less satisfactory 
in liquid Si at high temperature, where swift bonding and 
rcbonding processes take place. However, in all cases the 
dynamics is dissipative, most likely due to the fact that 
the propagation scheme employed is not symplcctic. It is 
possible to remedy this downward drift if we assume that 
the forces arising from our dynamics Fpc can be modeled 
as Fpc = Fbo — TflR/ which, as we shall see later, is an 
excellent assumption. The value of the intrinsic friction 
coefficient does not need to be known but it can be 
bootstrapped by taking a cue from the work of Krajewski 
and Parrinello 2l| . We sample the canonical distribution 



by using the following Langevin-type equation 
A/ / R / = F PC - 7i R / + H / , 



(7) 



where M/ is the ionic mass, jl is a Langevin friction 
coefficient and Bj = Bf + 3j a random noise. Using 
the assumption above this can be equally written as: 

M / R / = F bo -(7d+7l)R./ + 3/ (8) 

In order to guarantee an accurate sampling of the Boltz- 
mann distribution, the noise has to obey the fluctuation 
dissipation theorem: 

(3/ (0) B z (t)> = 6 { 1D + 7L ) MjkeTd (i) (9) 

The choice of is arbitrary, while the unknown has 
to be determined by requiring that the aggregate noise 
term generate the correct average temperature, as mea- 
sured by the equipartition theorem ^M/R^ = ^ksT. 
We will show that this leads to correct sampling. Since 
our initial dynamics is quite accurate is small, so that 
dynamical properties are also well reproduced. 

For the purpose of demonstrating our new method, 
we have implemented it in the mixed Gaussian Plane 
Wave (GPW) [22| code Quickstep j2j| which ispart of 
the publicly available suite of programs CP2K [24|. We 
present here calculations on metallic liquid silicon and 
liquid silica, to illustrate that our method works well irre- 
spective of band-gap, system size and system type. Both 
systems are known to be very difficult, and are exam- 
ples of liquid metals (Si) and of complex, highly polar- 
izable, ionic liquids (SiO^). In addition the simulations 
have been performed at 3000 K and 3500 K respectively, 
which leads to rapidly varying density matrix elements, 
thus making the propagation of the electronic degrees of 
freedom quite challenging. Hence, the selected tests can 
be considered as worst-case scenarios for our method. 

All simulations have been performed at their experi- 
mental liquid densities using double-zeta valence polar- 
ization (DZVP) basis sets, adequate densit y cu toffs, the 
Goedecker-Teter-Hutter pseudopotentials [25j and the 
local-density approximation to the exact exchange and 
correlation functional. For simplicity the Brillouin zone 
is sampled at the T-point only. Eq. ©is integrated using 
the algorithm of Ricci and Ciccotti [26| , with a time step 
of h = 1.0 fs. The friction coefficient 7l was set equal 
to zero, while the values for 7d turned out to be in the 
range of 10" 4 fs~ . We predict the new C's using K = 4 
in Eq. ([2]), which ensures time-reversibility up to Q(h 6 ). 
The OT minimizer is preconditioned, as proposed in [16j . 

We first consider the accuracy in terms of deviation 
from the BO surface. As can be seen in FIG. [1] our en- 
ergies are an upper bound to the ground state and are 
displaced by a small and approximately constant amount. 
It is also shown that the deviation from the BO surface 
can be even further reduced by increasing the number of 
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FIG. 1: Deviations from the BO surface of liquid S1O2 with 
respect to total energies (upper panel) and mean force devia- 
tions (lower panel). The deviation in the energies corresponds 
to a constant shift of 4.16 x 10 -4 Hartree per atom for one cor- 
rector step and 3.5 x 10 -5 Hartree per atom for two corrector 
steps. The average mean force deviation is unbiased. 
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FIG. 2: Partial pair-correlation functions g(r) of liquid Si 
(upper left panel) and liquid Si02 at 3000 K and 3500 K 
respectively, using a DZVP Gaussian basis set. 



corrector steps. In fact the deviation can be controlled by 
varying the number of corrector steps in order to achieve 
a preassigned accuracy level. In the following only simu- 
lations based on a single corrector step will be reported. 

We turn now to more realistic problems such as those 
shown in FIG. [2] Although these simulations have been 
performed with only a single corrector step, they are still 
amazingly close to the BOMD results. It should be em- 
phasized that even in liquid Si, which poses problems 
when using an ordinary Car-Parrinello scheme due to 
its metallicity, a single corrector step is sufficient. This 
establishes the efficiency of our method, since only a 
single preconditioned gradient calculation and no self- 
consistent iteration has to be performed. The possible 
acceleration, in comparison with regular BOMD calcula- 
tions, depends crucially on the system studied. In these 
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FIG. 3: The kinetic energy distribution calculated from a 
1 ns trajectory with a time step of 3.25 fs of metallic liquid 
Si64 using a DZVP Gaussian basis set and a density cutoff of 
100 Ry (left panel). Velocity autocorrelation function (upper 
right) and its Fourier transform (lower right) of 32 Water at 
325 K using a TZV2P Gaussian basis set, a density cutoff of 
280 Ry and the BLYP exchange-correlation functional. The 
Langevin friction coefficients are 7l = and ~ 10 -8 fs -1 . 



difficult cases a speed-up of two orders of magnitude com- 
pared to using the pure extrapolation scheme have been 
observed. For simpler problems an increase in efficiency 
of at least one order of magnitude can be expected. 

In FIG. [3] we present results which prove that also dy- 
namical properties can be evaluated with accuracy. This 
time we look at the velocity autocorrelation function and 
its Fourier transform at 325 K. The results are in good 
agreement with our reference calculations and are con- 
sistent with experiment and ab-initio all-electron calcula- 
tions (2?| , showing that in spite of the stochastic nature of 
Eq. JSJ) dynamical properties can also be simulated. This 
implies, that also non-equilibrium processes and chemi- 
cal reactions can be handled. In the same picture we 
verify that our assumptions are justified, and we are in- 
deed performing a canonical sampling, by showing that 
the kinetic energy distribution is Maxwellian. To this ef- 
fect, we have carried out a 64 atom liquid Si simulation 
for as long as 1 ns, to reduce the noise and to ensure a 
proper sampling of the kinetic energy distribution tails. 

Because of space considerations we have reported here 
only a fraction of the systems studied. In all cases our 
method has proven to be accurate and the gain in speed 
has always been remarkable. Also structure relaxation 
via annealing and geometry optimization have been suc- 
cessfully performed. Contrary to CPMD and related 
methods at least BOMD integration time steps can be 
used. Thanks to this development it is now possible to 
perform ab-initio molecular dynamics on medium-sized 
systems up to a few nanoseconds, thus making a new class 
of problems accessible to ab-initio simulations. The com- 
putational scaling of the algorithm is in principle linear. 
However, an efficient parallel sparse matrix algebra has 



not been implemented yet, so that the sustained scaling 
is still 0(MN 2 ), albeit with a much reduced prefactor. 

In conclusion we wish to highlight that our ideas are of 
general interest, as our method is independent of the par- 
ticular choice of the minimizer and propagation scheme, 
and can be further improved. In fact, during the sub- 
mission process we became aware of a time-reversible ex- 
trapolation scheme [28| , which could be used as an alter- 
native to ASPC. The possibility to apply the presented 
ideas with benefit to plane wave based CPMD is also to 
be underlined and will be presented elsewhere. 
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deep insights, and F. R. Krajewski and J. Kolafa for fruit- 
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and support from CSCS Manno, the ICT Services of ETH 
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